function hfile = vmr_LimitVMR(hfile, opts)
% VMR::LimitVMR  - apply intensity limits to V16 of VMR
%
% FORMAT:       [vmr = ] vmr.LimitVMR([opts])
%
% Input fields:
%
%       opts        1x1 struct with optional settings
%        .range     either 1x2 intensity range (limits included) or
%                   1x2 relative limiting (default: [0.001, 0.999])
%        .recalc8b  flag whether or not to recalc 8-bit portion {false}
%
% Output fields:
%
%       vmr         altered object
%
% Note: the VMRData16 field must be set correctly

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if isempty(bc.VMRData16) || ...
    numel(size(bc.VMRData16)) ~= numel(size(bc.VMRData)) || ...
    any(size(bc.VMRData16) ~= size(bc.VMRData))
    error( ...
        'BVQXfile:InvalidObject', ...
        'VMR16 data not found.' ...
    );
end
numvox = numel(bc.VMRData);
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'range') || ...
   ~isa(opts.range, 'double') || ...
    numel(opts.range) ~= 2 || ...
    any(isinf(opts.range) | isnan(opts.range) | opts.range < 0) || ...
    any(opts.range ~= fix(opts.range) & opts.range > 1)
    opts.range = [0.001, 0.999];
end
if ~isfield(opts, 'recalc8b') || ...
   ~islogical(opts.recalc8b) || ...
    numel(opts.recalc8b) ~= 1
    opts.recalc8b = false;
end

% interpret range
mmm = minmaxmean(bc.VMRData16);
vdt = single(bc.VMRData16);
mx = mmm(2);
opts.range = sort(opts.range(:));
if any(opts.range ~= fix(opts.range))
    hn = hist(vdt(:), single(0:mx));
    cn = cumsum(hn);
    if opts.range(1) ~= fix(opts.range(1))
        miny = find(cn > (numvox * opts.range(1)));
        if isempty(miny)
            opts.range(1) = 0;
        else
            opts.range(1) = miny(1) - 1;
        end
    end
    if opts.range(2) ~= fix(opts.range(2))
        maxy = find(cn < (numvox * opts.range(2)));
        if isempty(maxy)
            opts.range(2) = numel(hn);
        else
            opts.range(2) = maxy(end);
        end
    end
end

% limit VMR
minmask = [];
maxmask = [];
if opts.range(1) > 0
    minmask = (bc.VMRData16 < opts.range(1));
    bc.VMRData16(minmask) = opts.range(1);
    bc.VMRData16 = bc.VMRData16 - opts.range(1);
end
if opts.range(2) < mx
    maxmask = (bc.VMRData16 > opts.range(2));
    bc.VMRData16(maxmask) = opts.range(2) - opts.range(1);
end

% adapt VMR
if ~opts.recalc8b && ...
   (~isempty(minmask) || ...
    ~isempty(minmask))
    mmm = minmaxmean(bc.VMRData);
    if ~isempty(minmask)
        bc.VMRData(minmask) = mmm(1);
    end
    if ~isempty(maxmask)
        bc.VMRData(maxmask) = mmm(2);
    end
end

% set min/mean/max new
mmm = minmaxmean(bc.VMRData16);
bc.MinOriginalValue = mmm(1);
bc.MeanOriginalValue = round(mmm(3));
bc.MaxOriginalValue = mmm(2);

% try to perform good limitting
if opts.recalc8b
    vdt = single(bc.VMRData16);
    [hn, hx] = hist(vdt(:), single(0:bc.MaxOriginalValue));
    hx = uint16(hx);
    cn = cumsum(hn);
    lt = find(cn > 0.01 * cn(end));
    ut = find(cn > 0.99 * cn(end));
    if isempty(lt)
        lt = 1;
    else
        lt = lt(1);
    end
    if isempty(ut)
        ut = numel(cn);
    else
        ut = ut(1);
    end
    vdt(vdt > single(hx(ut))) = single(hx(ut));
    vdt = vdt - single(hx(lt));
    vdt(vdt < 0) = 0;
    mmm = minmaxmean(vdt);
    bc.VMRData = uint8(floor((225.99 / mmm(2)) .* vdt));
end

% set to content array
bvqxfile_setcont(hfile.L, bc);
